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ABSTRACT 

Scattering of light from an anisotropic source produces linear polarization in spectral lines and the 
continuum. In the outer layers of a stellar atmosphere the anisotropy of the radiation field is typically 
dominated by the radiation escaping away, but local horizontal fluctuations of the physical conditions 
£NJ , may also contribute, distorting the illumination and hence, the polarization pattern. Additionally, 

a magnetic field may perturb and modify the line scattering polarization signals through the Hanlc 
effect. Here, we study such symmetry-breaking effects. We develop a method to solve the transfer 
of polarized radiation in a scattering atmosphere with weak horizontal fluctuations of the opacity 
and source functions. It comprises linearization (small opacity fluctuations are assumed), reduction 
ly-j ■ to a quasi-planeparallel problem through harmonic analysis, and numerical solution by generalized 

standard techniques. We apply this method to study scattering polarization in atmospheres with 
horizontal fluctuations in the Planck function and opacity. We derive several very general results 
and constraints from considerations on the symmetries and dimensionality of the problem, and we 
give explicit solutions of a few illustrative problems of especial interest. For example, we show (a) 
. how the amplitudes of the fractional linear polarization signals change when considering increasingly 

smaller horizontal atmospheric inhomogeneities, (b) that in the presence of such inhomogeneities even 
a vertical magnetic field may modify the scattering line polarization, and (c) that forward scattering 
polarization may be produced without the need of an inclined magnetic field. These results are 
important to understand the physics of the problem and as benchmarks for multidimensional radiative 
' transfer codes. 

Subject headings: Polarization — Scattering — Line: formation — Sun: magnetic fields — Sun: 
atmosphere — Stars: atmospheres 
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> ' 1. INTRODUCTION 

oo ' 

The polarization generated by scattering in a spectral line depends sensitively on the global geometry of the scattering 
process (most notably, the distribution of incident radiation), as well as on the angular momentum of the atomic levels 
, involved in the transition. If the incident illumination is isotropic, the scattered radiation is unpolarized; if the 
• ■ illumination is collimated, the polarization of scattered light may be significant, reaching complete polarization at 90° 
' of the incident beam in a — > 1 transition (total angular momentum J = 0, 1 for the lower and upper levels of the 
, transition, respectively ). In the outer layers of a stellar atmosphere the radiation field is anisotropic mainly because 
light is escaping through the surface. When quantified (Sect. [5] below), the degree of anisotropy in a plane-parallel, 
semi-infinite atmosphere is typically of a few percent, which yields signals with polarization degrees on that order of 

> ! magnitude or lower (Stenflo et al. 1983a, b; Stenflo & Keller 1996, 1997; Gandorfer 2000, 2002, 2005). 

In the presence of a magnetic field the whole scattering process is perturbed and the ensuing polarization pattern 
altered. This phenomenon, the Hanle effect, is important because it can be exploited to detect and measure magnetic 
[ fields in astrophysical plasmas and, in particular, in the solar or a stellar atmosphere (e.g., Stenflo 1991; Trujillo Bueno 
- . i 2001; Casini & Landi Degl'Innocenti 2007). This goal requires a reliable modeling of the radiation field anisotropy that 
generates the polarization pattern in the first place. The anisotropy is dominated by the local center-to-limb variation 
of the radiation field, but horizontal fluctuations of the physical properties in the atmosphere may alter it too. 

In this paper we investigate the effect of horizontal fluctuations of the thermodynamic parameters on scattering 
polarization, how they compete with magnetic fields to modify the general polarization patterns of spectral lines, 
and how both effects may be disentangled for diagnostic purposes (see Manso Sainz & Trujillo Bueno 1999 for a first 
attempt at this program). We focus on weak opacity and Planck function fluctuations. This allows a semi-analytical 
treatment of the problem which greatly simplifies the numerical analysis and allows solving two and three dimensional 
problems with the computational cost of a plane-parallel one. 

We consider resonance scattering and the Hanle effect in a two-level atom whose lower-level is unpolarized; for 
definitcness, the results shown in this paper are for the Ji = — > J u = 1 transition. We also consider Rayleigh and 
Thompson scattering by exploiting its formal similarity to the resonance scattering problem. These problems can 
be treated within two alternative, but equivalent, formalisms. First, we may note that the polarization state of the 
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Fig. 1. — Left panel: reference system for polarization. A ray propagates along the direction CI, ei is perpendicular to f2 and lies on 
the meridian plane (the plane containing fi and the z-axis), and e2 is perpendicular to both, f2 and ei. Note that ei, e2, fi, form an 
orthonormal positive base. In Eqs. (J4j and JB), the positive Q-direction is taken to be along ei, i.e., perpendicular to the projected horizon 
(limb). 6g and \B are the polar and azimuthal angles of the magnetic field, respectively. Right panel: three lines-of-sight (LOS) are 
of particular relevance in our analysis for two-dimensional media (invariant direction along the y-eods). The slabs in the figure represent 
roughly the cosinusoidal fluctuations considered in Sect. 4. In case a, we observe close to the limb (fi = cos 9 = 0.1) and along the invariant 
direction (x = 90°); in case b, we observe close to the limb (fi = 0.1) and accross the slabs (\ = ); m case c, we observe at disk center 
( M = l). 

light irradiating the collectivity of atoms, and the polarization state of the scattered radiation are linearly related. 
Therefore, the whole scattering process may be described through a phase matrix relating the Stokes parameters of the 
incident and scattered beams; the source function of the radiative transfer equations is then an average over all possible 
incident angles of the phase matrix (e.g., Chandrasekhar 1960). Alternatively, we may consider the intermediate state 
involving atoms explicitly. The incident radiation excites the atoms, generating atomic level polarization (population 
imbalances and quantum coherences between sublevels); as a consequence, the reemited (scattered) radiation is linearly 
polarized. A complete description of the excitation state of the atomic system taking into account sublevel populations, 
coherences and statistical mixtures can be done through the density matrix (Fano 1957). Further, we will work with 
the spherical tensors components of the density matrix. In Sect. [5] we show how to formulate mathematically this 
problem using elements of this formalism. 

The density matrix formalism provides several advantages for the problem under consideration. First, its spherical 
tensor decomposition is especially well-suited to apprehend and exploit the symmetries of the problem. It will become 
apparent through the paper that symmetries lead to important simplifications in equations in a transparent way. 
Second, it is computationally advantageous to pose the problem with the atomic density matrix elements as unknowns, 
from which the radiation field straightforwardly derives, rather than to consider the radiation field itself as the unknown. 
This is because there are at most six density matrix components in our problem to be determined at each position 
in the medium, while there are three Stokes parameters that depend on the spatial, angular and frequency variables. 
Finally, the density matrix is a fundamental and very general concept and polarization transfer theories based on it 
have been developed to treat more general problems than the ones considered here (Landi Degl'Innocenti & Landolfi 
2004), e.g., resonance scattering with lower- level atomic polarization, zero- field dichroism (Trujillo Bueno & Landi 
Degl'Innocenti 1997; Trujillo Bueno 1999; Trujillo Bueno et al. 2002; Manso Sainz & Trujillo Bueno 2003), and 
scattering in multilevel systems (Manso Sainz & Landi Degl'Innocenti 2002; Manso Sainz, Landi Degl'Innocenti, & 
Trujillo Bueno 2006; Manso Sainz & Trujillo Bueno 2003, 2010; Stepan & Trujillo Bueno 2010, 2011). Therefore, the 
methods and techniques developed here can then be generalized to these other problems. 

Coherent continuum scattering polarization, resonance scattering polarization, and the Hanle effect in plane-parallel 
media have been studied extensively. On the contrary, in two and three dimensions, the problem has been attacked more 
sparingly (but see Dittmann 1999; Paletou, Bommier, & Faurobert-Scholl 1999; Manso Sainz & Trujillo Bueno 1999; 
Manso Sainz 2002; Trujillo Bueno, Shchukina & Asensio Ramos 2004; Trujillo Bueno & Schukina 2007, 2009; Schukina 
& Trujillo Bueno 2011; Anusha & Nagendra 2011; Anusha, Nagendra, & Paletou 2011). This is because the inherent 
computational complexity of multidimensional radiative transfer problems increases notably when polarization is taken 
into account, due to its sensitivity to the radiation field angular dependence. Here, we overcome these difficulties by 
considering scattering polarization transfer in atmospheres that are only weakly inhomogeneous, in the sense that 
horizontal fluctuations can be considered as a perturbation of a plane-parallel problem. We linearize and derive 
equations for the fluctuations of the radiation field (Sect. [3]). Then, harmonic analysis of the linear problem is possible 
(Sect. |4|). This approach allows application of simple symmetry arguments to simplify the problems and isolation of 
the relevant variables in each case, which leads to a system of transfer problems that are formally similar to the plane- 
parallel case and can thus be solved applying standard numerical techniques (Sect. [5]). This program generalizes the 
approach of Kneer & Heasley (1979), Kneer (1981), and Trujillo Bueno & Kneer (1990) for the unpolarized problem. 

2. FORMULATION AND GENERAL EQUATIONS 

This section states mathematically three polarization transfer problems: resonance scattering (Sect. 12. ip . coherent 
continuum scattering (Sect. 12. 2j) . and the Hanle effect (Sect. I2.3j) in optically thick atmospheres. General symmetries 
of these problems mainly related to their dimensionality are naturally described within the formalism, and discussed 
in Sect. El 



Scattering polarization in inhomogeneous atmospheres 



3 



2.1. Scattering line polarization 

We consider resonance scattering polarization in an inhomogeneous, non-magnetized, static atmosphere. In such 
a medium scattering produces linear polarization but not circular polarization. Since the atmosphere is horizontally 
inhomogeneous, all three Stokes parameters /, Q and U are necessary to describe the polarization state of the light. 
The lower-level of the transition is assumed to be unpolarized (all its magnetic sublevels are equally populated). 
Therefore, there is no selective absorption of polarization components (i.e., no zero- held dichroism) as discussed in 
Trujillo Bueno & Landi Degl'Innocenti (1997), Trujillo Bueno (1999) and Manso Sainz & Trujillo Bueno (2003). The 
transfer equations for polarization then read: 

±I = - K (I-S T ), (la) 
±Q = - K (Q-S Q ), (lb) 
d U = -k(U - Su), (lc) 



d 



s 



where d/ds is the derivative along the ray path, k is the absorption coefficient, and Si, Sq and Su the source functions 
for the corresponding Stokes parameters. In the presence of a background continuum the absorption coefficient has 
contributions from the line and the continuum: k — k e <p v + n cont (</>„ is the line absorption profile), while the source 
function 5/ = r I/ S , J ine + (l-r v )S'J ont , with r v = k^^/k. We will assume the continuum to be in LTE, S c j ont = B V {T), 
(B V (T) is the Planck function) and unpolarized (Sq — r v S l Q nc , Su = r^S 1 ™ ). A polarized continuum due to Rayleigh 
or Thompson scattering can be readily included without difficulty (see next subsection). 

The line source functions depend on the excitation state of the atoms which, within the density matrix formalism, 
is described through the density matrix p (Fano 1957; Blum 1981). Its diagonal elements p m m are the populations of 
the individual magnetic sublevels and the non-diagonal ones the coherences between them. It is convenient to express 
the density matrix in its spherical components pq (K = 0, 2 J, Q = —K, K; Brink & Satchler 1968). Introducing 

the new variables Sq — ^pr~ ^2J + +i PQ ( v * s ^ e f rec l uenc y 01 the transition and Mt the total population of the 
ground level) we will derive expressions that are a natural generalization of the standard radiative transfer theory for 
unpolarized radiation. The Sq components are, in general, complex quantities satisfying the conjugation property 
[Sq]* = (—1)^S_q ("*" stands for complex conjugation), as a consequence of the hermiticity of the density matrix p. 
For computational purposes it is more convenient to work with the following real quantities 



% = n[S§ + i-l) Q S? Q ], (2) 



$q=^[Sq -(-i) Q s5 Q L (3) 

defined for Q > (the equivalent Q < elements are redundant). Note that Sq and Sq arc the real and imaginary 
parts of Sq* (Q >Q). 

The line source functions can be obtained applying the density matrix theory of polarization in spectral lines (e.g., 
Landi Degl'Innocenti & Landolfi 2004) assuming complete frequency redistribution. They can be expressed as 

Sf lc = S° + w <£ I-L(3f - l)S 2 - VSp^T^^osxSf - sinxSf) 

V3. ^ ^ 



^-(1 - M 2 )(cos2 X 5| - sin2 X 5 2 2 )} 
) f 3 

n 
2 



(4b) 

X (l + M 2 )(cos2 X 5 2 2 -sin2x4 2 )}, 



= Vs{ v / W 2 (sin X ^ 1 2 + cosx^r 2 ) 

. . , (4c) 

+ ^(sin2x^| + cos2x^ 2 2 )). 

where 9 = arccos p and x are the inclination and azimuth of the ray, respectively, and the reference directions for 
polarization are shown in Fig. Q] (i.e., the positive-Q direction is along ei). The coefficient Wj J( is a numerical factor 
that only depends on the quantum numbers of the transition considered (see Eq. (10.12) in Landi Degl'Innocenti & 
Landolfi 2004). Explicit values for several transitions, as well as the limits for large J, are compiled in Table 1. 
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Equations Q express the angle- dependent line source functions in terms of the six (angle and frequency independent) 
variables S$, S§, Sf, Sf, Sf, S 2 satisfying (e.g., Trujillo Bueno & Manso Sainz 1999): 

(5a) 
(5b) 
(5c) 
(5d) 
(5e) 
(5f) 

In Eqs. (J5J) , e = C u z/{C u t + A u g) is the collisional destruction probability {C u i is the collisional desexcitation rate 
and A U £ the Einstein coefficient for spontaneous emission), and 5^ = /A u i is the depolarizing rate due to elastic 
collisions with neutral hydrogen normalized to A u g. Finally, the radiation held tensor components read: 
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(6a) 
(6b) 
(6c) 
(6d) 
(6e) 
(6f) 



where df2' = sin6*'d#'dx' = d/i'dx'. The Jq and Jq components with Q — 1,2, in Eqs. (JSJ) correspond to the real and 
imaginary parts, respectively, of the complex spherical component Jq of the radiation tensors discussed in Sect. 5.1 of 
Landi Degl'Innocenti & Landolfi (2004). Assuming complete frequency redistribution, the actual Jq components to 
be introduced in Eqs. (J5J) are the frequency averages over the absorption profile: 



J$ = I d^JgM- 



(7) 



Hereafter, we will suppress the subscripts expressing explicitly the dependence of the different variables on v, unless 
necessary to avoid ambiguities. 

In the limit 5^ -> 00, the S Q elements vanish (see Eqs. p^ -(|5f])). Then S% ne = S[} ne = 0, S^ 110 = (Eqs. 0}), 
and we recover the well-known expression for a two level atom under complete frequency redistribution (Eq. I5a[) of the 
classical radiative transfer theory for unpolarized radiation (e.g., Mihalas 1978). 

We note in passing, that the equivalent phase matrix formalism is recovered by writing explicitly the expressions for 
the radiation field tensors (Eqs. ©-©) into Eqs. ([5]), and then into the expressions for the source functions (Eqs. (@|). 
Thus proceeding we obtain explicit expressions for the source functions of the Stokes parameters, at a given angle and 
frequency, as averages over the incident radiation field. The linear operator relating incident and emergent radiation 
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is the so-called phase matrix (e.g., Chandrasekhar 1960). The methods developed here are thus straightforwardly 
extended to the phase matrix formalism. 

2.2. Coherent continuum scattering polarization 

Scattering by free electrons (Thompson scattering) and by atoms (most notably, neutral hydrogen) in the far wings 
of resonance lines (Rayleigh scattering) is coherent in the electron/atomic rest frame (Dirac 1925), but otherwise, 
the scattering phase matrix is equivalent to that of resonance scattering in a Ji = — > J u = 1 transition (e.g., 

Chandrasekhar 1960). Therefore, we may write w^q = 1 in Eqs. (J3J), e = 6^ = in Eqs. (JS|), and consider Eqs. as 
a convenient factorization of the phase matrix in terms of six variables Jq that are angle independent. 

For passing to the laboratory frame, we must consider the Doppler effect of the Maxwellian distribution of the 
scatterers velocities. This implies convolving the incident radiation field, at each frequency v 1 and angle //, with a 
profile 

<!>(&■ v,v') = ; 1 exp(- , ^"' / j, t , 1, (8) 

V ' A D ^(l-cose) I 2(1 - cos 6) A 2 , J w 

where 6 = arccos(-y/l — /i 2 \/l — n' 2 cos(x — \') + MA*') i s the angle between incident and scattered beams, and Ad 
the Doppler width corresponding to the scatterers (note, in passing, that Au(e~) 3> Ad(H°)). 

If the incident continuum radiation field is spectrally flat, a convolution over v' with the expression in Eq. ([5]) leaves 
the spectrum unaltered. Scattering can thus be considered coherent also in the laboratory frame, and continuum 

(2) 

scattering is described by exactly the same equations of the previous subsection with the formal substitutions w 10 = 1 , 
e = 5^ = 0, <p v = 8{y — v'). Such an approximation for the continuum is formally equivalent to the complete frequency 
redistribution approximation for line scattering by atoms. 

It is important to recall that this approximation neglects important effects arising from the presence of spectral 
structure in the incident radiation field. In particular, such an approximation should be considered with care (if at 
all) when interpreting the continuum in the vicinity of spectral lines (see Landi Degl'Innocenti & Landolfi 2004). 

The source functions 5'| catt , Sq att and Sfj att due to pure coherent Rayleigh or Thompson scattering in the continuum 

are given by Eqs. ©-([I]) with the following formal substitutions: Wj 2 j ( = 1, e = <5^ 2 ^ = 0. Let a be the absorption 
coefficient for scattering, K cont the true continuum absorption coefficient due to thermal absorption, k = a + K cont , and 
s = a/K, then the total source functions are Si = sS s j catt + (1 - s)B v , S Q = sS S Q att , and S v = sSff att . 

2.3. Hanle effect 

The Hanle effect is the modification of scattering line polarization due to the presence of weak magnetic fields. By 
weak magnetic field we mean the regime in which the splitting of the emission and absorption profiles (~ z/l, the 
Larmor frequency), is of the order of the natural width of the upper level of the transition (i.e., A U £ 4 ), and hence, 
negligible compared to the Doppler line width Ar/D- In this limit, the splitting of the line profile in its a and w 
components is negligible and the only noticeable effect of the magnetic field on the polarization state of light stems 
from its perturbing the excitation and coherence state of the atoms. Therefore, the radiative transfer equations are still 
given by Eqs. ([T]) and (|4]). However, the excitation state of the atomic system is now a competition between radiative 
excitation and partial decoherence induced by the magnetic field splitting (besides collisional excitation/relaxation). 
The statistical equilibrium equations for the Sq components read (e.g., Manso Sainz & Trujillo Bueno 1999) 
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In Eqs. (|9b[) - ([9fj) . T = 0.879 x 10 7 g u B/A u e (g u is the upper-level Lande factor and B the magnetic field strength in 
gauss), s a = sin a, c a = cos a, and 8b and xb are the inclination and azimuth of the magnetic field, respectively (see 

The micro-turbulent limit of these equations is often studied because it leads to simpler equations and it is considered 
a suitable approximation to interpret spectropolarimetric observations with low spatio-temporal resolution (Stenflo 

4 Note that we are assuming that the lower level is unpolarizcd and the so-called lower-level Hanle effect regime is not considered in this 
paper. 



G 




Fig. 2. — A two-dimensional scattering medium is invariant under translations along a direction (say, the y-axis), and reflections on the 
x-z plane. Thus, at an arbitrary point in the atmosphere the polarization of light propagating along two rays fl and f2' is obtained by 
reflection from the other. Light polarized along directions ei and in ray f2, is polarized along directions and e^, respectively, in ray 

fi'. Hence, If_i :X = ^m>-x anc ^ Qm,x = — x (Eqs. Illlal l- llllbl l). However, light polarized along the direction (ei ± e2)/\/2 in ray fl is 
polarized along direction (e^ =p e' 2 )/\/2 in ray O' and U^^ x = — £/ Mi _ x (Eq. micf ). Note that, therefore, U = at right angles from the 
fluctuations (x = 0) and, in particular, at disk center observation (fi = 1). 

1982, 1994; Landi Degl'Innocenti 1985). In this limit, one assumes that the magnetic field changes its orientation 
and strength at scales much smaller than the photon mean free path. The excitation state of the atom at a point in 
the atmosphere is then the average over all possible realizations of the magnetic field and the Sq components satisfy 
averaged Eqs. © (once formally solved). Averaging over an isotropic distribution of fields of constant strength gives 
equations identical to those of the pure scattering case (Eqs. ([5])), but with the Jq radiation field tensor components 
multiplied by a simple numerical factor (Trujillo Bueno & Manso Sainz 1999) 

where 7 = T(l - e)/[l + 6^(1 - e)}. 

Note that the same symmetry considerations given for scattering polarization in Sect. [2"7T1 hold in the presence of an 
isotropic micro-turbulent field. 



2.4. Symmetry considerations 

Equations ([TJ, together with Eqs. (HJ-© are general and describe the transfer of polarized line radiation in a 
scattering atmosphere of arbitrary geometry. However, the presence of symmetries in the medium leads to important 
simplifications. Thus, in a plane-parallel atmosphere, symmetry imposes that light be linearly polarized either parallel 
to the horizon or perpendicular to it. Choosing these directions as the reference for polarization (see Fig. (Jl])), the 
only non- vanishing Stokes parameters are / and Q, and they are independent of the azimuthal angle. Consequently, 
Jq = for Q ^ 0, which implies Sq^q — 0, and we recover the well-known expressions for scattering line polarization 

in a plane-parallel medium (e.g., Trujillo Bueno & Manso Sainz 1999). 

Consider a general Cartesian two-dimensional medium. We choose without loss of generality the y-axis as the 
invariant direction. Then, the following symmetries hold (see Fig. [2] for proof): 

I/j,,x — Ip,-xi (H a ) 

Qn, x = Q^.-xi ( llb ) 

U^x = -U^-x- (He) 

Introducing these expressions into Eqs. (j6|), one finds that J\ = Jf = 0. Therefore, S\ = if = 0, and just four 

variables (Sq, Sq, Sf, and Sf) are necessary to completely describe the excitation state of the atomic system in a 
two-dimensional medium. 

The symmetries of the problem can be further exploited, but this requires a more explicit knowledge of the hori- 
zontal fluctuations of the physical system (Sect.[4j. Clearly, the same symmetry considerations apply to the coherent 
continuum scattering polarization case as well. 

In general, the presence of a deterministic magnetic field breaks the symmetry of the radiation field. Thus, regardless 
of the problem's dimensionality three Stokes parameters (/, Q, and U) are necessary to describe the polarization state 
of the radiation field, which is no longer axially symmetric. Therefore, the Jq elements do not vanish in general and 
all six Sq elements are necessary to describe the excitation state of the atomic system. The description simplifies only 
for very special geometrical configurations. Thus, if the medium is planeparallel and the magnetic field vertical, we 
recover the symmetries of the scattering polarization in a planeparallel medium discussed in Sect. 12.11 In fact, that 
problem is completely equivalent to the scattering polarization problem: the Hanle effect does not operate. 
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In a two-dimensional medium with no magnetic field, just four Jq (Sq ) components are necessary to describe the 
illumination (excitation state) in the medium (Sect. [XT] Manso Sainz 2002). This holds true also in the presence of a 
horizontal magnetic field aligned along the invariant direction since in that case the symmetries illustrated in Fig. [2] 
are also satisfied. Note, however, that unlike the planeparallel case, a vertical magnetic field does modify scattering 
polarization and hence the Hanle effect does not vanish. 

3. WEAK FLUCTUATIONS AND LINEARIZATION 

3.1. Radiative transfer equations 

Now, we consider a weakly inhomogeneous atmosphere along the horizontal direction. We assume that all the 
physical variables can be expressed as a small perturbation of a plane-parallel atmosphere and that we can thus 
neglect second order terms. We write 

k(x, y, z) = fi(z) + 8k(x, y, z), (12) 
Si(x,y,z) = Si(z) + SSi(x,y,z), (13) 
I{x,y,z) = I[z) + 5I(x,y,z), (14) 

(with analogous expressions for Q and U), where / is the solution of the planeparallel transfer problem posed by R 
and Si. (To avoid lengthy and heavy algebraic expressions, throughout the paper we will drop explicit dependencies 
on the independent variables — space, frequency, angles — whenever they are not strictly necessary and does not lead 
to ambiguities.) Introducing expressions (fT2)) - (IT4l) into the radiative transfer equations (UJ we get 

—51 = -k(6I - 5Si) - 5k(I + 61 - Si - 5SA, (15a) 
as 

^-5Q = -R(5Q - SS Q ) - 5k{Q + 6Q-S q - 5S q ), (15b) 

^-5U = -k(5U - 5S V ) - 5k(U + 5U - S v - 5S V ). (15c) 
as 

If 6k = 0, Eqs. (|15|) are linear, irrespectively of the magnitude of the amplitude of the source function fluctuation. In 
the presence of small opacity perturbations (5k/ k -C 1), and assuming that the ensuing perturbation on the intensity 
and linear polarization is small (6I/I,6Q/I,6U/I <C 1), we neglect second order terms in Eqs. (fT"5|) : 

-^51= -R(6I-6Sf), (16a) 
as 

.^-6Q = -k(6Q-6S c q s ), (16b) 

A<5£7 = -k(6U - 5Sf), (16c) 
as 

5Sf = SSj-— (I -Si), (17a) 

6S$ - 6S Q - S 4(Q-S Q ), (17b) 

SSff = 6Su- S 4(U-Su). (17c) 

K 

Equations t|l 6[) - f| 1 T|) are linear in the sense that if 5±I and are solutions of transfer problems with 5\K, 5±Si, and 
62^, 52S1, respectively, then aSiI + $82! is the solution of the transfer problem with fluctuations ab~\K + j352K and 
aSiSi + /3S2S1, and analogously for Q and U. Note in passing, that due to linearity, a problem with fluctuations on 
both, opacity and source function, can be retrieved from the superposition of a problem in which only the opacity 
fluctuates (while the source function remains constant), and the problem in which only the source function fluctuates 
(while the opacity is fixed). 

We will consider the following boundary conditions. An open boundary with no incident illumination at the top 
(z = Zm)' 5I^\(x,y, zm) = 8Q fl \(x,y 1 zm) — 8U^\(x,y, Zm) — 0, for /1 < 0; a thermal, unpolarized radiation field at 
the bottom (z — z m ), well below the thermalization depth: 5I^\(x, y, z m ) = 5B V , 5Q fl \(x. y, z m ) ~ 5U^\(x, y, z m ) = 
for fi > 0. 

3.2. Statistical equilibrium equations 

The statistical equilibrium equations for all the problems considered in the previous section can be expressed as an 
algebraic system of equations: 

MS = b, (18) 



where the effective source functions are 



where S is a formal vector whose (up to 6) elements are the Sq components, M is a square matrix whose coefficients 
depend on the parameters of the problem (magnetic field strength and orientation, collisional rates e, 5), and the 
components of the formal vector b are linear combinations of Jq components and the Planck function. 
Here we consider fluctuations in the Planck function and opacity: 

B V = B V + 5B V , (19) 
k = R + 6k. (20) 

This implies fluctuations for Sq and hence, for the radiation field tensors, which are decomposed as 

S$ = S$ + 5S%, (21) 

J$ = J$ + 6J$. (22) 

We thus find that if the mean unperturbed variables satisfy the statistical equilibrium equations MS — b, the pertur- 
bations satisfy the same system M8S = db. 

It is possible to consider more general fluctuations in other physical parameters of the problem, as for example, in 
the magnetic field or in the collisional rates. It is then possible to linearize the statistical equilibrium equations which 
yields the following formal equations for the perturbations 

M8S = Sb SMS. (23) 

4. HARMONIC ANALYSIS 

4.1. Radiative transfer equation 

We perform an harmonic analysis of the linear problem posed in the previous section and we show that the general 
multidimensional radiative transfer problem posed by Eqs. (|16p ~ (|17p reduces to several coupled plane-parallel radiative 
transfer problems. We consider sinusoidal fluctuations along the x-axis of the opacity k and Planck function B u of the 
form 

B(x,y,z) = B(z) + AB(z)coskx, (24) 
k(x, y, z) = k(z)[1 + acos kx], (25) 

If a > 0, the fluctuations are said to be in phase; if a < 0, the fluctuations are said to be in anti-phase. For the 
linearization considerations of the previous section to apply, we assume \a\ <C 1. 
Due to the linearity of the problem the Stokes parameters have the following spatial dependence 5 

SI (x, y, z) = Ail(z) cos kx + A 2 I(z) sin kx, (26a) 

SQ(x, y, z) = A\Q(z) cos kx + A 2 Q{z) sin kx, (26b) 

5U(x, y, z) = AiU(z) cos kx + A2U(z) sin kx, (26c) 

and 

8Sf(x,y,z)^[A 1 Si(z) - a(I - Si)]coskx + A 2 Si(z)smkx, (27a) 

6S Q s {x,y,z) = [AiSq(z) - a(Q - Sq)] cos kx + A 2 S Q (z) sin kx, (27b) 

6Sf{x,y,z) = [AxSuiz) - a(U - S v )]coskx + A 2 Su(z) sinkx. (27c) 

The amplitudes Ai i2 S7, Ai^Sq, and Ai^SV, are expressed in terms of Ai^Sq, Ai^Sq,..., Ai^Sf , as in Eqs. (U]). 

Now, we derive radiative transfer equations for the amplitudes A1.2/, Ai^Q, and Ai^U as follows. In Cartesian 
geometry, taking the y-axis as the invariant direction, the explicit expression for the d/ds operator in the radiative 
transfer equation is 

d _ d + x d ^ 
ds dz dx ' 

with /i = cos 6 and A = sin 6 cos x (see Fig.H}. From Eqs. (fT6]) - (fi"?]l and Eqs. (J26J) - ((2HJ) we get, after some elementary 
algebra, the following coupled system of radiative transfer equations for the amplitudes: 

^A,I = A,I - [AiS, - a(I - S z ) - ^A 2 7] , (29) 
dr L k 

d \k 
li—A 2 I = A 2 I - [A 2 Si + —Ail] , (30) 
dr k 

5 Higher harmonics cos(nfca:) and sin(nfca:), with n > 2, appear only if they are already present on the radiation field illuminating the 
boundaries. But then, lacking sources, they are exponentially attenuated and in a semi-infinite mean atmosphere, they vanish completely 
when approaching the surface. 
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Fig. 3. — The atmosphere defined by Eqs. 124M251 is invariant under a 180° rotation around the origin of coordinates. Thus, radiation 
along a ray f2 at point P(x, y, z), is identical to radiation along ray fl' at point P(—x, —y, z). Therefore, Ifj, x ( x ) = ^x+» ( — x )> Q^x( x ) = 
Qfj, x +Tr{-x), and U^ x (x) = x). 

with analogous expressions for Ai i2 Q and Ai. 2 [/. A change of variable has been introduced in Eqs. (j2"9")) - (l3T)l) . from the 
geometrical z-scale to the the monochromatic optical depth scale r, which is dehncd from the mean opacity k, such 
that dr = — Rdz. Note that unlike a true plane-parallel problem, here, the radiation field depends explicitly on the 
azimuthal angle. The apparent asymmetry between Eq. (|29p and Eq. pOj) (a term on a appears in the former but not 
on the latter), is just due to our choice of cosinusoidal opacity fluctuations in Eq. (|25j) . If a more general dependence 
were chosen (e.g., cos(fcir + </>)), an analogous term would then appear in Eq. (f30|) too. 

The boundary conditions for the whole fluctuations imply to the following boundary conditions for the amplitudes: 
Ai^x^m) = &i,2Qhx( z m) = &i,2U mx (zm) = 0, for (j, < 0, at the top; Ai/ MX (z m ) = AB(z n 
Ai^Q^ x (z m ) = Ai ! 2^ Atx (z m ) = 0, for fj, > 0, at the bottom. 



-nth A 2 -/^ x (z m ) 



4.2. Statistical equilibrium equations 
From Eqs. ()26[) . the spatial dependence of the radiation field tensor perturbations is 

SJq = Ai Jq (z) cosfcx + A 2 Jq (z) sinfcx, 



(31) 



where the amplitudes Ai. 2 J^f (z) are computed from the corresponding Ai. 2 J, Ai^Q, and Ai^U amplitudes according 
to Eqs. (j6|). Equation (|3ip . in turn, implies the same spatial factorization for the SSq fluctuations, 

5Sq = A 1 S%(z)coskx + A 2 S%{z)smkx, (32) 

as can be seen from the statistical equilibrium equations of the problem under consideration. Moreover, due to the 
linearity of the problem, the very same statistical equilibrium equations are satisfied, independently, for the terms 
fluctuating as coskx (A\Jq, A\Sq), and the terms fluctuating as sinfcx (A 2 Jq^ A%Sq). 

4.3. Further symmetry considerations 

It is possible to exploit the additional symmetries arising in the problem from the particular fluctuation law we 
are considering (Eqs. (f2"4 |) - (f25"|) ). From Fig. [3l it is clear that SI l _ lx (x) = 5Ip X +„(—z), SQ^x) — <5Q AIX+W (— x), and 
8U^ x {x) = <5f7 Alx + 7r (— x), which implies 

AiI MX = AiI MX+ ,r, (33) 
A 2 / MX - -A 2 / A1X+T , (34) 

(the same is true for A12Q and Ai j2 [7). Equations (|33|) - (l34l) lead to the following straightforward identities: 



dx A 2 / x = / d X Ai/ X sinx = / dx A 2 / x sin2x = 0, 
Jo Jo 



(35) 



which apply to A\^Q and Ai 2 C7 as well. When these relationships are taken into account in the definitions of the 
radiation field tensors (Eqs. ([6|), we get 



<5J ° = AiJ °(z)cos(fca;), 
SJq = Ai Jq(z) cos(kx), 
<5J 2 = A 2 J 1 2 (z)sin(fcx), 



<5J 2 



Ai J 2 (z) cos(kx). 



(36a) 
(36b) 
(36c) 
(36d) 
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This factorization of the radiation held tensors is a consequence of the symmetry illustrated in Fig. @. This symmetry 
does not hold anymore in the presence of an inclined magnetic field; the radiation field tensors then fluctuate according 
to the more general Eqs. (j3"Tj) . Note, in passing, that this would also be the case, even without magnetic fields, if the 
perturbed variables had a phase difference between them (i.e., other than or n), e.g., SB ~ cosfcx and 6k ~ sinfcx. 

The most important point to keep in mind is that the general multi-dimensional problem has been reduced to a set 
of plane-parallel coupled ones given by Eqs. (|29|) -(|30 j) and their corresponding set of statistical equilibrium equations. 
These quasi-planeparallel problems must now be solved numerically. 

5. NUMERICAL METHOD OF SOLUTION 

We iterate to obtain the self-consistent solution of the radiative transfer equations for the Stokes parameters ampli- 
tudes and the statistical equilibrium equations for the density matrix amplitudes. Section [5.11 shows two different, but 
related, methods to integrate the RT equations to obtain the Stokes parameters amplitudes and the radiation field 
tensor amplitudes A 12 </q when the Ai 2 Sq elements are known. Section 15.21 in turn, shows how to calculate new 

values of the A\ }2 Sq amplitudes from the just computed radiation field amplitudes. We explain how to iterate this 
process to reach convergence. 

5.1. Formal solution of the transfer equation 
5.1.1. Method (a) 

We integrate the RT equations along short characteristics (SC). This method is based on the following strategy. Let 
M, O, and P be three consecutive points along a ray path. Assuming that the source function varies parabolically 
between them, the intensity at point O may be expressed as 

7(0) = I(M)e- A ™ + <p M Si(M) + ipoSriO) + ^Sj(P), (37) 

where Atm is the optical distance between M and O, and ipM, 4>o and if>p are coefficients that only depend on 
the optical distance between points O and M (Atm), and M and P (Arp). The explicit functional form of the ip 
coefficients for linear or parabolic variations of the source function between grid points can be found in Kunasz & 
Auer (1988). The linear version (without the ij)p term) must be used at the boundaries when no P point is available. 
Expressions analogous to Eq. (1571) are valid for the Stokes parameters Q and U. In planeparallel problems, M, O, 
and P correspond to three consecutive points on the z-grid. In two- and three-dimensional problems, M and P do 
not correspond, in general, to grid points and some interpolation is needed (see Kunasz & Auer 1988; Auer & Paletou 
1994; Auer, Fabiani Bendicho, & Trujillo Bueno 1994; Fabiani Bendicho & Trujillo Bueno 1999). Here, we avoid 
this computationally intensive interpolations by exploiting the known functional dependence of the variables on the 
horizontal direction. 

Let [xltVl^ z l) be the Cartesian coordinates of point L. Noting that AiI(zl) = 51(0, Ul, z l): and A 2 I(zl) = 
5I($-,y L ,z L ), then Eqs. QKft-QTQ, and Eq. fl37J imply 6 



2k 

Afc Xk 
A 1 I(z ) = [Ai/(z A f)cos( — A M ) - A 2 /(z M )sin( — A M ) 

fi fi 



-Atm 



Xk Xk 
+ ^m[A 1 Si(z m )cos(—A m ) ~ A 2 S I (z M )sm( — A M )] +^ A 1 Si(z ) (38) 
A 4 M 

Xk Xk 
+ 4> P [AiSV(zp) cos(— A P ) + A 2 S I (z P ) sin(— A P )1 , 
H M 

Xk Xk 
A 2 I(zo) = [AJIzm) sin(— A m ) + A 2 J(z M ) cos( — A M )]e~ A ™ 

fi fi 

Xk Xk 
+ ^ M [A 1 S I {z M )sm( — A M ) + A 2 S7(z M ) cos( — A M )] +^ A 2 Si(z ) (39) 
1 fi fi 

Xk Xk 
+ tpp \ - AxSi(zp) sin( — A P ) + A 2 Si(z P ) cos( — A P )1 , 
fi fi 

and analogously for Ai 2 Q and Ai^ 2 U, with Am = zo — Zm and Ap — zp — zo- Note that the Stokes parameters 
amplitudes A\ >2 I, Ai^Q, and A\ >2 U depend on just one spatial dimension (depth), but they depend on both angular 
variables, fi — cos 9 and A = sin 9 cos x- Proceeding from one boundary to the other along the ray direction, Eqs. (|38)) - 
(|39|) allow us to calculate A\ f2 I (and Ai i2 Q, and Ai t2 U) at each height of the vertical grid, and from them the Stokes 
parameters at any other horizontal position if necessary. 

5.1.2. Method (b) 

Equations (l2"9"l) - (|30p are formally identical to the RT equation in a plane-parallel atmosphere if we identify the 
term between brackets on the right hand side with a source function. We exploit this formal similarity to integrate 

6 For the sake of clarity, only the case without opacity fluctuations (a = 0) is treated here explicitly. The formal substitution A\Sj — > 
AiS — a(I — Si) in these equations yields the general expressions with opacity fluctuations. 
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numerically these equations through the well-known short-characteristics (SC) method for the formal solution of the 
RT equation. Thus, integrating explicitly Eqs. (j29|) - (l30|) between three consecutive spatial points M, O and P along 
the ray assuming a parabolic variation for the terms between brackets, the intensity amplitudes at point O can be 
expressed as 6 

Xk 

AiI(*o) - A 1 /(z M )e- A ™ +^ M [A 1 S I (z M ) - — -A 2 I(z M )] 

K[ZM > (40) 
AA; Xk 
+ Vo[AiSj(z ) - —-^-A 2 I(z )] +ipp[AiSi(z P ) - — — A 2 /(z P )], 

K{Zo) K(Zp) 

Xk 

A 2 I(z ) = A 2 I(z M )e- A ™ +^ M [A 2 Si(z M ) + tA x 1{z m )\ 

<zm) , } 

Afc AA: 

+ Vo \A 2 Sl(zo) + rrr-rA^izo)] + V>P \A 2 Si(z P ) + ___AiI(*p)] , 
K[Zo) K{Zp) 



with similar expressions for Ai^ 2 Q and Ai^U. Note that here, as in Eqs. (|38|) - (I39|) . Atm is calculated from the mean 
opacity R. 

In contrast with a true plane-parallel problem, Eqs. (|4T))) - (l4"T1) do not give explicitly the values of A X I and A 2 I at a 
given point O once the intensity at the previous point M is known. Instead, Ai/ and A 2 I at O depend on those very 
same quantities at points M and P. The intensity amplitudes are given implicitly, in the form of an algebraic system 
of equations that can be written in matrix form as 

- A v M + B v - C v P = b Q , (42) 



where 



Co - 



0j 




b _ ( iPmAiSi(zm) +4>oAiSi(z ) + tppAiSi(zp) 
° y iPmA 2 Si(z m ) + ^oA 2 Si{z ) + ippA 2 Si(zp) 

and vm, fo and vp are vectors whose two components are the amplitudes A\I and A 2 I at points M, O and P, 
respectively. By changing the corresponding values in vectors b and v, the expressions apply to all Stokes parameters. 
The set of Eqs. (j42|) for all points in the atmosphere form a 2x2-blocks tridiagonal system of equations. Such a system 
can be efficiently solved through a standard forth and back substitution scheme (e.g., Press et al. 1992): 

Vi = DiV i+1 + u h (43) 

Di = (Bi- A i D i _ 1 )~ 1 C i , (44) 

Ui = (Bi - AiDi-x)- 1 ^ + Aiiti-t), (45) 

where i = 1,...,N Z (N z is the number of spatial grid points). We start at one atmosphere's boundary (i = 1) and 
calculate Di and Ui at all grid points until i = N z — 1. At the other boundary i = N z , Cn z = Dm, =0 and 
vn, = UNz ■ Then we proceed backwards obtaining Vi from equation (|43p . This process involves the inversion of 2 x 2 
matrices which can be performed analytically, and only requires a little extra numerical effort. 

Both methods of integration of the RT equation presented in this section converge to the same solution for increasingly 
finer grids (i.e., as Az — ¥ 0). Method (a) assumes a parabolic variation of the source functions perturbations 8Si(x, y, z) 
between three consecutive points, regardless of the known factorization given in Eqs. (|26p . Method (b) assumes a 
parabolic variation of the source functions amplitudes A\^ 2 Si{z) between three consecutive heights, and the horizontal 
variation of the physical parameters is taken into account explicitly on the derivation of the equations. Therefore, 
method (a) is easier to implement and slightly faster (a tridiagonal system of equations is avoided), while method 
(b) is more accurate for coarser spatial grids. The difference between both is more evident for higher k, and very 
inclined rays close to the y-z plane. In that case, rays may cross several horizontal periods of the fluctuating source 
functions, which is poorly approximated by a parabolic variation. Therefore, method (a) requires finer z-grids at high 
wavenumbers k. 

5.2. Iterative scheme 

Equations (|3"5)) - (f3T))) or Eqs. (|4T)l) - (|4Tj) show a linear relationship between the source functions and Stokes parameters 
amplitudes. Formally: 

81= A n [5Sj] + 5X, (46a) 
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5J° 


= Aoo^Sq - 




f- A 02 5S 2 - 


h Am6§* - 


h A 0c <5S COIlt - 


\-8JS 


8J% 


= A 10 SS° - 


h A n 5Sg - 


FAi2<5S?- 


h A 14 <5S 2 - 


h A lc 5S cont - 




SJl 


= A 20 SS° - 


h A 21 5S* - 


H A 22 <5S 2 - 


h A 24 <5S 2 - 






8Jl 


= AaaSS* - 


h - 










8Jl 


= A iQ SS° - 


h A 4 i5Sg - 


H A 42 5S 2 - 


h A 44 <5£ 2 " 






SJl 


= A 53 S§1 - 


h Ag 5 5S| - 











SQ= A n [SS Q ]+SQ, (46b) 
5Z7= A n [5S u ]+SU ) (46c) 

where <5J, and 8U are 2 x N z - element vectors whose components are the Ai i2 I, Ai j2 <5, Ai. 2 t/ amplitudes at 
the A z heights in the atmosphere, respectively. The 2xN z components ol 51, SQ, and 8U are the corresponding 
Stokes parameters amplitudes transmitted from the boundary. SSi, SSq, and SSjj are 2xiV z -elenient vectors with 
components Ai i2 SV, Ai j2 Sq, and Ai^Sjj, respectively. Finally, Ajj is a (2 x N z ) x (2 x 7V 2 )-element matrix; the 
same for all three Stokes parameters that depends on both, the inclination (through /i) and azimuthal angle (through 
A) of the ray. The An(i,j) elements are involved functions of the optical distances between spatial grid points that 
need not be evaluated explicitly — they are implicitly calculated every time a formal solution is performed as in the 
previous section. However, Eqs. (|46|) show that they can be numerically evaluated straightforwardly. If the medium 
is not illuminated at the boundaries {SI = 8Q = 8U = 0), and we consider a source function perturbation vanishing 
at all points but at point i, where SSi(i) — 1 (or SSg(i) — 1 or SSu{i) — 1), then the corresponding Stokes parameter 
fluctuation calculated at point j as in Sect. 15. H is 5l(j) = An(i,j) {8Q(j) = An(i,j), 8U(j) = An(i, j)). In particular, 
the diagonal element A x ti(i,i) is the intensity amplitude at a given point when the corresponding source function 
amplitude at that point is 1 and vanish elsewhere. 
We introduce the 2 x A z -element vectors SSq, SSq, SS^, and SJ°, SJq, SJ% , whose components are Ai i2 Sq, 
Ai^S"! , Ai j2 Jq, Ai j2 </|, respectively, at each point in the atmosphere. The source functions depend linearly on 
the Sq components (Eqs. Q), while the Jq elements are weighted angular averages of the Stokes parameters (Eqs.|HJ). 
Formally: 

(47a) 
(47b) 
(47c) 
(47d) 
(47e) 
(47f) 

where the matrices A a p are obtained from Afj: 

Kp(i,j)= Jdv4>u(i)r v (j) y^wjfAo(i,j), (48) 

Aoc(i,j) - / di/&,(i)[l - r v (j)]I^^ c A n (i,j). (49) 

The angular factors w^f and are obtained after some simple algebraic manipulations (see Appendix) and arc 
given in Tables [5] and [6j respectively. Note that due to symmetry arguments raised previously, Eqs. (|47dl) and (|47fl) 
are relevant only if a magnetic field is present. 

Let SSq° be the value of SSq at some previous "old" iterative step and let SJq° be the vector of radiation field 
tensor amplitudes calculated from SSq ° W by formal integration of the RT equation (Sect. 15. 1[) . We seek to calculate 
the "new" values SSq " cw of the density matrix amplitudes at the following iterative step. We proceed by a variation 
of the approximate A-iteration (ALI) method along the lines used by Manso Sainz & Trujillo Bueno (1999), and Manso 
Sainz (2002). The new amplitude value at a given point "i" in the atmosphere is obtained by solving the statistical 
equilibrium equations of the problem (Eqs. ©or ([5])) with the following estimates for the radiation field tensors: 

« <5J °M° ld + Aoo(*, *)[^ °(z) n ™ - SS^) ol % (50) 
8J 2 Q {i)^8Jl^r d . (51) 
Substituting Eq. ((SDJ into Eq. fl5aj (or ([9"a])t leads to 

scO, -view _ (1 ~ e)[^o°(») old - A 00 (i,i)JSg(»)° ld ] +eB„ 

6S ° {1) ~ l-(l- e )Aoo(i,i) ' (5) 

which is the well-known ALI correction (see Trujillo Bueno & Manso Sainz 1999, and references therein). Substituting 

Eq. (HU into Eqs. (|5b))-([5f| or Eqs. d9bj) ) - dHJ) leads to a A-iteration correction for SSfi, SSf, S 2 . 

This is the simplest convergent iterative scheme possible for this problem. The iterative scheme can be generalized 
by extending the procedure on Aoo in Eq. ([50|) to any other operator A a/ g in Eq. [37] This leads to slightly more 
complex expressions of the type of Eq. (1521 for all SSq elements. However, there is little gain in the convergence rate 
for the particular problems we considered here. 
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Fig. 4. — Maximum relative change Rq(Sq) (left panel) and Rc(Sq) (right panel), for different values of the horizontal wavenumber k 
(labels). Scattering polarization in a non magnetized atmosphere with constant properties (and identical for all cases) is considered. Dotted 
line: convergence rate for k = with Ng-acceleration. 



Figure I5.2I shows the maximum relative change 



Rn = max 



ssg(i) now -asg(i) old 

SS K (j)„ow I 



(53) 



of Sq and Sq in a scattering polarization problem. Different horizontal wavenumbers are considered. It can be shown 
that the convergence rate improves for larger horizontal wavenumbers. This is because effectively, the medium becomes 
more and more optically thin. It can also be seen that a factor ^3 improvement is obtained by just performing Ng- 
acceleration to the standard ALI iterative scheme. Comparing both panels we see that the convergence rate of the Sq 
elements is tied to the convergence of the populations (Sq). 
The convergence rate is independent of the method (a or b) employed for the formal solution integration. 



5.3. Numerical Considerations 
5.3.1. Angular quadrature and symmetry 

The angular quadrature in multidimensional radiative transfer problems requires special consideration. This is even 
more so when dealing with polarization. 

Kneer and Heasley (1979) introduced an efficient numerical quadrature to deal with the complicated angular de- 
pendencies of the radiation field in the unpolarized limit of the problems considered here (sinusoidal fluctuations of 
opacity and Planck function) . It consists of a Gaussian quadrature for polar angles taken from the invariant direction 
(here, the y-axis), and a trapezoidal rule for the azimuth on the orthogonal plane. Unfortunately, this quadrature is 
no well suited for polarization purposes. In order to capture the symmetries of the problem and to be able to recover 
the plane-parallel limit (both, when k — and k — > oo), the azimuthal quadrature must satisfy the following relations 
exactly: 

E ; cos \i = S-i sin Xi = cos 2xi = sin 2xi = 0, (54) 

where the sum extends over all the directions with a fixed inclination with respect to the vertical. When Eqs. (|54j) are 
not satisfied exactly, spurious sources of polarization appear. 

The simplest quadrature satisfying these properties is a trapezoidal rule with at least four points uniformly distributed 
between and 2tt. We use a iV M -point Gaussian quadrature for the inclination and a ./V x -point trapezoidal rule for the 
azimuth. As pointed out by Kneer and Heasley (1979), this choice has the inconvenience of not resolving adequately 
the variations of the radiation field due to sinusoidal fluctuations of the source function. Consequently, relatively fine 
grids are needed at moderate k to avoid oscillations very close to the surface, where the medium becomes optically 
thin. This inconvenience is overcome by our not introducing spurious sources of polarization. 

The precision of the finite grid discretization is quantified by 

/a ir \ F [Av] - F[ Sv ] 

e F (Av\6v) = j (55) 



which gives the value of some quantity F in a coarse grid of a given variable (Au) relative to its value in a finer grid 
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Fig. 5.— Relative discretization error for the value of AJ§/B„ (O), AJ^/J^ (A), Ajf/jfi (□), and Ajf/J§ (x), at z = (solid lines), 
and z = 2 (dotted lines), for an atmosphere with horizontal fluctuations with k = 0.3 of the source function. The relative errors are given 
with respect to the solution given in Table 3. Panel a: frequency discretization, N v is the number of frequency points between the core 
and 10 Ai/q Panel b: polar angle discretization, Ng is the number of Gaussian points Panel c: azimuthal angle discretization, N x is the 
number of points for a trapezoidal rule between and 2n. Panel d: spatial discretization, N&z is the number of grid points between z and 
2+1. 

(Sv), all other quantities kept equal. Panels b and d of Fig. 15.31 show the eF(A9\S9) and €f(Ax\Sx) for F — AJq/B v , 
AJ 2 /J °, AJf/Jo, and AJf /J °, at z = and z = 2. The reference grid is the one for the results given in Fig. 8 and 
Table 3. 



5.3.2. Frequency quadrature 

We consider Gaussian line profiles <j) v and a static medium. A trapezoidal rule is used for the frequency integration. 
Panel a in Fig. 15.31 show ep{Av\8v) as a function of the number of points N v between v\ — (line core), and 
vm u = lOAi/p (Ai/rj being the Doppler width). 

5.3.3. Spatial (vertical) discretization 

In plane-parallel media only the optical depth is a relevant variable, the actual distribution of opacity with height 
is unimportant; in two and three dimensions however, the stratification of opacity with geometrical depth becomes 
relevant too. The atmosphere's gravitational stratification suggests an exponential mean opacity R(z) = e~ z . Here, 
K is normalized to its value at z — and the spatial variables z and x are given in units of the opacity scale height 
H. For reference, in the solar photosphere H is between w 100 and 200 km for typical spectral lines. The (vertical) 
optical depth (see Sect. I4.ip is then r = e~ z — e ~ 2max . 

We consider a spatial grid {zk}k=i,...,N z between z m - m = —12 and z max = 18. Therefore, r s» e~ z in units of H; 
hence, r min = 10~ 8 and r max = 1.6 x 10 5 . Panel c in Fig. l5.3l show the relative error ep(Az\Sz) with respect to the grid 
necessary to obtain the results in Fig. 8 and Table 3 as a function of the number of z grid points per unit z interval. 

5.3.4. Accuracy and performance 

When we consider horizontal fluctuations only on the Planck function the radiative transfer problem posed is already 
linear and the methods of solution presented in this section (in particular, the formal solution) are exact. Furthermore, 
since no horizontal discretization is needed, the solution obtained is more accurate than the solution obtained with 
other methods (based either on finite differences or SC) that require the discretization of all three spatial coordinates 
(e.g., Auer, Fabiani Bcndicho, & Trujillo Bucno 1994; Manso Sainz & Trujillo Bucno 1999). Thus, for example, such 
'classical' methods approach the solutions of the problems considered in Sects. 6.1 and 6.2 below only in the limit of 
very fine horizontal grids. 
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Treating the horizontal fluctuations exactly has the additional advantage of improving the numerical performance 
by avoiding costly interpolations. SC methods require knowledge of the source function at points M and P, and the 
intensity at point M (see Eq. (|37[) . But for very special configurations, in two- and three-dimensional spatial grids, the 
M and P points do not belong to the spatial grid and the relevant quantities must be interpolated from neighboring 
grid point values. This requires that the interpolation coefficients be calculated for every direction of the angular 
quadrature. This calculation scales at best 7 as ~ Na x Np (iVo = x N x being the total number of directions of 
the angular quadrature; Np — N x , or N x x N y in two and three dimensions, respectively, N x and N y being the total 
number of points along the respective axes) , and should be repeated at each iterative step (everytime a formal solution 
is performed), or otherwise, stored in the computer's memory. 

When the opacity fluctuates horizontally, the radiative transfer problem is linearized along the lines of Sect. 3. 
The numerical considerations just discussed apply to the linearized problem too: accuracy and performance improve 
due to the lack of horizontal discretization. However, there is an intrinsic loss of precission in the solution from the 
linearization of RT equation. The relative error in the solution is of the order of the second order terms neglected in 
the linearization: ad I for the amplitude SI (and analogously for the other Stokes parameters amplitudes). Typical 
errors in the amplitudes of the profiles for the problems considered in Sects. 6.3 and 6.4 are below 1-0.1%. 

6. RESULTS AND DISCUSSION 

We study a few selected, idealized problems: coherent continuum scattering polarization in the presence of horizontal 
fluctuations of the temperature (Sect. 16. 1| . resonance scattering polarization with horizontal fluctuations of the 
temperature (Sect. 16. 2|) . resonance scattering polarization with horizontal fluctuations of the opacity and source 
function (Sect. 16. 3|) . and Hanle effect with horizontal fluctuations of the opacity and source function (Sect. I6.4j) . In 
particular, we consider cosinusoidal fluctuations of the parameters and we study the behavior of the solutions and 
emergent Stokes polarization profiles with the wavenumber. 

6.1. Continuum scattering polarization, source function fluctuations 

We consider coherent scattering polarization in a constant properties, semi-infinite atmosphere with sinusoidal, 
horizontal fluctuations of the Planck function. 

The total opacity (absorption plus scattering) is exponentially stratified with height. We normalize spatial distances 
to the opacity scale height such that k = k c + a — n e~ z , with n cont /cr constant throughout the atmosphere. An 
arbitrary height reference is set by choosing k = 1. The optical depth from the surface (z — > oo), at a given fi is 
then r = e~ z /[j,; conversely, the height (in opacity scale units) at which r = 1 along a line-of-sight (LOS) with /i is 
z = log(l//zv / 7r). Finally, we normalize AB, AJq, A5q amplitudes to the mean Planck function. 

The solution of this problem for K cont / (K cont + a) = 10 -4 and a Planck function fluctuation AB/B = 1 is given in 
Fig. [6] and Table [2] They show the variation of AJq, AJq, AJf, and AJ| with height as a function of the horizontal 
wavenumber k. Due to linearity, the solution for any value of AB/B is obtained as the linear combination of the 
plane-parallel case (here, the k = solution) and AB/B times this one. 

For symmetry reasons (see Sect. 12.41) . 5 J\ = 6J$ — and therefore, SSf = <5<Sf = 0; moreover, Ai Jf = A 2 J| = 
(see Sect. EM and therefore, A ± Sf = A 2 5| = 0. 

The behavior of AJq and A5q is only slightly affected by polarization and it closely resembles the behavior of the 
mean intensity and source function fluctuations of the unpolarized radiation problem (cf. Kneer 1981). 

On the other hand, the amplitudes AJq, AJf, and AJ| (hence, the corresponding A5q, ASf, and ASf ones), 
share the following asymptotic behaviors with depth (z — > — oo), and close to the surface (z — ¥ oo). Deeper than 

the thermalization depth, the radiation field becomes isotropic and therefore, AJq,AJ^,AJ| — > 0; towards the 
surface, horizontal transfer smears out horizontal fluctuations, thus recovering the unperturbed limit and therefore, 
AJq,AJ^,AJ| — > (but for k — 0, which corresponds to the planeparallel limit). Between these two limits, the 
amplitudes rise in absolute value, reach a maximum at some height, and decrease again. This maximum absolute 
value tends to zero as k — ¥ oo. The reason is that for increasingly small structures, they uniformly become optically 
thin, approaching the planeparallel unperturbed limit {AJq — > 0) throughout the whole atmosphere. 

It is interesting to study the behavior with wavenumber k at a given height, for example, at z ~ 2.3, which 
corresponds to r„ rs 1 at /i = 0.1 (i.e., close to the limb). The amplitudes monotonically increase (in absolute value) 
with wavenumber up to k ~ 0.2; beyond that, AJq and ASq become smaller the larger the horizontal wavenumber k 
— which will explain the behavior of the emergent fractional polarization (see below). When horizontal inhomogeneities 
appear, they break the rotational symmetry of the radiation field, which perturb the (vertical) anisotropy AJq and 
yields Jf , Jf 7^ 0. Then, at some point, those horizontal fluctuations become optically thin at the 'height of formation' 
and finally, the radiation field recovers its original symmetry. 

Figure [7] shows the center-to-limb variation of the Stokes profiles amplitudes when observing along the 'slabs', 
i.e., along the j/-axis. Due to symmetry reasons, for observations along the invariant direction, / and Q fluctuate 
cosinusoidally, (i.e., A 2 I = A 2 Q = 0) just like the perturbation, while U fluctuates sinusoidally (i.e., AiU = 0). This 

7 For very inclined rays, or when the horizontal grid is relatively fine compared to the vertical grid, long-characteristics integrations 
might be necessary, worsening the problem of interpolation by so doing across several different cells for every single ray (e.g., Auer et al. 
1994). 
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Fig. 6. — Polarization by coherent scattering in an atmosphere with horizontal, cosinusoidal fluctuations (labels indicate wavenumber k) 
of the Planck function (Eq. H24H *), An exponentially stratified atmosphere with ft cont /(<7 + K cont ) = 10 -4 , AB/B = 1 is considered. From 
top to bottom, and from left to right: AJ^, AJq/JJ, AJ^/jO, AJ^/Jq, respectively, as a function of height (z) measured in units of the 
opacity scale height (z = corresponds to r = 1). 

can also be understood by taking x = 90° (A = 0) in Eqs. (1^1) - ([501) . Then, amplitudes Ai and A 2 are independent 
between them, and the just mentioned result follows from Eqs. evaluated at x = 90°, and the fact that the only 
non vanishing components are A5°, ASq, ASf, and A5|- 

The limb darkening of the A\I amplitudes adds up to the limb darkening of the mean model. At large k, however, 
the fluctuations become optically thin, AI — > at all fi, and the limb darkening law of the unperturbed model is 
recovered. Amplitudes AiQ and A2U increase towards the limb before steeply falling to zero at /i — > 0; no horizontal 
transfer effects exist in the limit of tangential observation. This is a consequence of the amplitudes ASq falling towards 
zero when approaching the surface. The non monotonic behavior of ASq with varying wavenumber k at a fixed height 
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TABLE 2 

Radiation field tensor amplitudes for K cont /(K cont + a) = 10~ 4 



AJO/B 



k 



z 


T 


0.1 


0.2 


0.3 


0.5 


0.7 


1.0 


8 


0.00033 


5.50(-3) 


2.36(-3) 


9.1(-4) 


l.l(-4) 


1.3(-5) 




6 


0.0025 


6.59(-3) 


3.44(-3) 


1.62(-3) 


3.1(-4) 


5.3(-5) 




1 


0.018 


8.12(-3) 


5.13(-3) 


2.94(-3) 


8.3(-4) 


2.2(-4) 


3.0(-5) 


2 


0.135 


1.16(-2) 


8.69(-3) 


5.99(-3) 


2.48(-3) 


9.5(-4) 


2.3(-4) 





1.0 


2.71(-2) 


2.31(-2) 


1.85(-2) 


1.07(-2) 


5.81(-3) 


2.31(-3) 


-2 


7.4 


1.271-1) 


1.18-1) 


1.06(-1) 


7.91(-2) 


5.67(-2) 


3.37(-2) 


-4 


54.6 


6.12(-1) 


5.98(-l) 


5.78(-l) 


5.26(-l) 


4.69(-l) 


3.88(-l) 


-6 


403.4 


9.99(-l) 


9.98(-l) 


9.97(-l) 


9.93(-l) 


9.88(-l) 


9.77(-l) 





k 



z 


T 


0.1 


0.2 


0.3 


0.5 


0.7 


1.0 


8 
6 
1 
2 

-2 


0.00033 

0.0025 

0.018 

0.135 

1.0 

7.4 


1.61(-1) 
1.50(-1) 
1.241-1) 
7.14(-2) 
9.32(-3) 
-5.4(-4) 


1.27(-1) 
1.39(-1) 
1.33(-1) 
8.66(-2) 
1.53(-2) 
-3.2(-4) 


7.26(-2) 
9.65(-2) 
l.ll(-l) 
8.70(-2) 
2.02(-2) 
-l.l(-5) 


1.54(-2) 
3.06(-2) 
5.25(-2) 
6.14(-2) 
2.33(-2) 
6.2(-4) 


2.52(-3) 
7.49(-3) 
1.93(-2) 
3.38(-2) 
1.99(-2) 
1.08(-3) 


1.4(-4) 

7.9(-4) 

3.65(-3) 

1.16(-2) 

1.24(-2) 

1.35(-3) 


AJ2/J0 


z 


T 


0.1 


0.2 


k 

0.3 


0.5 


0.7 


1.0 


8 
6 
1 
2 

-2 


0.00033 

0.0025 

0.018 

0.135 

1.0 

7.4 


-1.09(-1) 
-9.55(-2) 
-7.20(-2) 
-3.851-2) 
-8.38(-3) 
-2.6(-4) 


-9.63(-2) 
-1.03(-1) 
-9.50(-2) 
-6.11(-2) 
-1.54(-2) 
-5.3(-4) 


-5.68(-2) 
-7.46(-2) 
-8.43(-2) 
-6.62(-2) 
-1.99(-2) 
-7.9(-4) 


-1.23(-2) 
-2.43(-2) 
-4.14(-2) 
-4.89(-2) 
-2.16(-2) 
-1.23(-3) 


-2.02(-3) 
-6.02(-3) 
-1.54(-2) 
-2.73(-2) 
-1.80(-2) 
-1.49(-3) 


-1.2(-4) 

-6.4(-4) 

-2.97(-3) 

-9.47(-3) 

-1.10(-2) 

-1.57(-3) 




z 


T 


0.1 


0.2 


k 

0.3 


0.5 


0.7 


1.0 


8 
6 
4 
2 

-2 


0.00033 

0.0025 

0.018 

0.135 

1.0 

7.4 


-4.39(-2) 
-3.60(-2) 
-2.34(-2) 
-9.27(-3) 
-1.59(-3) 
-3.8(-5) 


-4.23(-2) 
-4.351-2) 
-3.70(-2) 
-2.03(-2) 
-4.56(-3) 
-1.4(-4) 


-2.58(-2) 
-3.30(-2) 
-3.54(-2) 
-2.51(-2) 
-7.15(-3) 
-2.7(-4) 


-5.75(-3) 
-1.12(-2) 
-1.86(-2) 
-2.09(-2) 
-9.20(-3) 
-5.5(-4) 


-9.7(-4) 

-2.85(-3) 

-7.19(-3) 

-1.24(-2) 

-8.22(-3) 

-7.5(-4) 


-5.4(-5) 

-3.1(-4) 

-1.46(-3) 

-4.59(-3) 

-5.39(-3) 

-8.6(-4) 



manifests here as an increase and decrease of Q and U with A: at a given \i. 

6.2. Scattering line polarization with source function fluctuations 

The main difference between coherent and resonance scattering is that in the latter case, photons may redistribute 
(here, completely) within the spectral profile. As a consequence, the thermalization depth is shallower in the former 
case (e.g., Mihalas 1978), which leads to steeper gradients of the source function and hence, to larger radiation field 
anisotropies and emergent fractional polarization. As it turns out, horizontal transfer effects — and hence, A Jq — are 
also reduced in the resonance case. 

We consider resonance scattering in a constant properties, semi-infinite atmosphere with sinusoidal, horizontal 
fluctuations of the Planck function. In particular, we consider an infinitely strong line («; cont = 0). The line opacity 
exponentially stratified with height. We normalize spatial scales to the opacity scale height (k 11110 = n^ nc e~ z ), and an 
arbitrary height reference is set by choosing k 1 " 1c = 1. We consider a Gaussian absorption profile (thermal motions 
and microturbulence dominate), and we normalize frequencies to the thermal Doppler width (v = (y— v®)/ /Avd] 4>v = 

e~ u /\/tt)- The optical depth at the line center, taken from the surface (z —> oo), at a given [i is then r = e~ z /^/ttii] 
conversely, the height (in units of the opacity scale) at which r = 1 along a LOS with /j, is z = log(l / ' ■ Finally, 
we normalize AB, AJq , ASq amplitudes to the mean Planck function. 

Figure [5] and Table |3] show the solution of this transfer problem for a strong scattering line (K cont = 0, e — 10 -4 ), 
and a Planck function fluctuation AB/B = 1. They show the variation of AJq, AJq, AJf, and AJf, and the 

corresponding AS®, ASq, ASf, and ASf, with height as a function of the horizontal wavenumber k. Clearly, the 
linearity and symmetry considerations given in Sect. (|6.1[) apply here too. 
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Fig. 7. — Center-to-limb variations of amplitudes A\I/I, AiQ/I, and A2U/I corresponding to the amplitudes in Fig. (??), for observation 
centered and along the slabs, i.e., at x = 0, and along the y axis (LOS with \ = 0)- Labels indicate wavcnumber k. 

As in the coherent case, A J® and ASq are only slightly affected by polarization if compared with the solution of 
the unpolarized problem (6 — > 00). As seen by the k = limit in Figs. ?? and [3 the 2-level atom thermalizes deeper 
in the atmosphere (Mihalas 1978). The shallower gradient of the source function leads to a smaller anisotropy (Jgj 
e.g., Manso Sainz 2002), and to smaller horizontal anisotropies (AJq). Otherwise, the behavior of A Jq with height is 
identical in the coherent and CRD case. 

It is interesting to study the behavior with wavenumber A: at a given height, say, z w 1.73, which corresponds 
to t„ w 1 at /1 = 0.1 (i.e., close to the limb). The amplitudes increase (in absolute value) monotonically with 
wavenumber up to k ~ 0.2; beyond that, AJq and ASq become smaller the larger the horizontal wavenumber k - 
which will explain the behavior of the emergent profiles seen in Fig. (fT0|) (see below). First, horizontal inhomogeneities 

perturb the (vertical) anisotropy AJq, breaking the rotational symmetry of the radiation field, which yields Jf , Jf ^ 0. 
Then, at some point, horizontal fluctuations become optically thin at the 'height of formation' and finally, the radiation 
field recovers its original symmetry. 

From the statistical tensors amplitudes ASq we compute the amplitudes of the emergent polarization profiles for 
any direction. In particular, Fig. [9] shows the center-to-limb variation of the line core Stokes profiles amplitudes 
when observing along the 'slabs', i.e., along the y-axis. As in Sect. 16.11 the behavior with k of the emergent Stokes 
parameters amplitudes at a given /i can be understood from the behavior with wavenumber of ASq at a given height 
(remembering that, for example, z sa 1.73 corresponds to r„ w 1 at fi = 0.1). 

Due to symmetry reasons, for observations along the invariant direction, / and Q fluctuate cosinusoidally, (i.e., 
A 2 I = A 2 Q = 0) just like the perturbation, while U fluctuates sinusoidally (i.e., AiU = 0). This can also be 
understood by taking \ = 90° (A = 0) in Eqs. (|2"9"|) -([30 |) . Then, amplitudes Ai and A2 are independent between them, 
and the just mentioned result follows from Eqs. ((4]) evaluated at \ = 90°, and the fact that the only non vanishing 
components are AS$, AS$, ASf, and ASf . 

Figure [TO] shows the full profile of the emergent fractional polarization for a specific fluctuation with AB / B = 
0.2. This has been calculated using the linearity of the problem, from the amplitudes of the Stokes profiles of the 
planeparallel case (k = 0), and 0.2 times the amplitude of the sinusoidally perturbed model. When observing at 
fj, = 0.1 and along the y-axis (left panel), the amount of polarization increases with wavenumber for k < 0.2; at k ~ 0.2 
it changes trend, and for larger k values the plane-parallel limit is again recovered. A similar behavior is observed at 
disk center (right panel) , where the presence of horizontal fluctuations alone breakes the symmetry of the problem and 
is able to generate a polarization signal. Since there is more radiation along the slabs than perpendicularly to them, 
the polarization direction is perpendicular to the slab (i.e., along the x-axis). 



6.3. Scattering line polarization, source function and opacity fluctuations 

We consider now the additional effect of an opacity fluctuation. If the opacity and Planck function are in phase, 
hotter regions are more opaque, cooler regions more transparent, and radiation is channeled from brighter to fainter 
regions (Cannon 1970). If, on the contrary, opacity and Planck function fluctuate in anti-phase, hot regions are more 
transparent (hence, brighter), emission in cooler regions is further blocked (thus fainting), and contrast is enhanced. 

The effect of these mechanisms on the emergent polarization is illustrated in Fig. [TT] It shows the effect of opacity 
fluctuations in phase (a = +0.1), or in anti-phase (a = —0.1) with the Planck function, for fluctuations with k = 0.3, 
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Fig. 8. — Scattering polarization in an atmosphere with horizontal, cosinusoidal fluctuations of the Planck function (Eq. 1241 ). Labels 
indicate wavenumber k. A — > 1 transition with e = 10 — 4 and no depolarizing collisions is considered. From top to bottom, and from left 
to right: Ajg, Ajg/Jg, AJ2/J0, AJ 2 2 /J°, respectively (solid lines), and ASg, ASg/Sg, AS^/Sg, ASf/Sg (dashed lines). 

and AB/B = 0.2 (as in Fig. [lOj The modification of the Q/I profiles by such opacity inhomogeneities is negligible 
(~0.1%), compared with the intrinsic polarization from the unperturbed model and the Planck function fluctuation 
itself in close-to-the-limb observations (fj, — 0.1). More interesting is the polarization signal at disk center (/i = 1), 
which is exclusively due to the horizontal inhomogeneities. Opacity fluctuations with a > (a < 0) enhance (diminish) 
the polarization signal generated by the Planck function fluctuations. 

Close to the limb, anti-phase fluctuations (a < 0) hamper blurring of horizontal inhomogeneities, which enhances 
Stokes U at the edge between hot and cool slabs (x = -P/4). Fluctuations in phase enhance radiation channeling, 
hence horizontal blurring, and Stokes U decreases. At disk center (fi = 1), U = for symmetry reasons (see Fig. 2). 
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TABLE 3 

Radiation field tensor amplitudes for e = 10 -4 



AJO/B 



z 


r 


0.1 


0.2 


0.3 


0.5 


0.7 


1.0 


8 


0.00033 


4.21(-3) 


1.30(-3) 


3.64(-4) 


2.62(-5) 


1.5(-6) 




6 


0.0025 


5.07(-3) 


1.91(-3) 


6.55(-4) 


7.24(-5) 


8.4(-6) 


1.9(-7) 


1 


0.018 


6.181-3) 


2.84(-3) 


1.19(-3) 


1.99(-4) 


3.65(-5) 


4.7(-6) 


2 


0.135 


8.01(-3) 


4.45(-3) 


2.27(-3) 


5.68(-4) 


1.57(-4) 


3.51(-5) 





1.0 


1.33(-2) 


8.72(-3) 


5.34(-3) 


1.94(-3) 


7.70(-4) 


2.63(-4) 


-2 


7.4 


3.501-2) 


2.62(-2) 


1.861-2) 


9.22(-3) 


4.87(-3) 


2.27(-3) 


-4 


54.6 


1.13(-1) 


9.43(-2) 


7.61 (-2) 


4.87(-2) 


3.221-2) 


1.92(-2) 


-6 


403.4 


3.33(-l) 


3.04(-l) 


2.71(-1) 


2.14(-1) 


1.70(-1) 


1.25(-1) 


Ajg/Jg 


z 


r 


0.1 


0.2 


k 

0.3 


0.5 


0.7 


1.0 



0.00033 

0.0025 

0.018 

0.135 

1.0 

7.4 



1.41(-1) 
1.36(-l) 
1.22(-1) 
8.89(-2) 
3.29(-2) 
1.92(-3) 



8.33(-2) 
9.66(-2) 
1.03(-1) 
8.90(-2) 
4.29(-2) 
5.87(-3) 



3.50(-2) 
4.94(-2) 
6.38(-2) 
6.74(-2) 
4.08(-2) 
7.87(-3) 



4.34(-3) 
9.17(-3) 
1.77(-2) 
2.78(-2) 
2.53(-2) 
7.67(^3) 



4.95(-4) 
1.55(-3) 
4.43(-3) 
1.02(-2) 
1.34(-2) 
5.76(-3) 



3.9(-5) 

1.4(-4 

6.6(-4) 

2.56(-3) 

5.43(-3) 

3.46(-3) 



z 


r 


0.1 


0.2 


0.3 




0.5 


0.7 


1.0 


8 


0.00033 


-1.04(-1) 


-6.52(-2) 


-2.78(- 


-2) 


-3.47(-3) 


-3.8(-4) 


-l.l(-5) 


6 


0.0025 


-9.72(-2) 


-7.48(-2) 


-3.91(- 


-2) 


-7.36(-3) 


-1.24(-3) 


-l.l(-4) 


4 


0.018 


-8.331-2) 


-7.85(-2) 


-5.03(- 


-2) 


-1.42(-2) 


-3.60(-3) 


-5.5(-4) 


2 


0.135 


-5.95(-2) 


-6.81(-2) 


-5.33(- 


-2) 


-2.25(-2) 


-8.38(-3) 


-2.15(-3) 





1.0 


-2.80(-2) 


-3.79(-2) 


-3.56(- 


-2) 


-2.18(-2) 


-1.16(-2) 


-4.75(-3) 


-2 


7.4 


-5.54(-3) 


-8.58(-3) 


-9.43(- 


-3) 


-8.00(-3) 


-5.76(-3) 


-3.41(-3) 


AJI/JO 


z 


r 


0.1 


0.2 


0.3 


k 


0.5 


0.7 


1.0 



0.00033 

0.0025 

0.018 

0.135 

1.0 

7.4 



-4.33(-2) 
-3.89(-2) 
-3.07(-2) 
-1.87(-2) 
-7.30(-3) 
-1.44(-3) 



-2.92(-2) 
-3.27(-2) 
-3.29(-2) 
-2.62(-2) 
-1.32(-2) 
-3.09(-3) 



-1.28(-2) 
-1.77(-2) 
-2.22(-2) 
-2.22(-2) 
-1.41(-2) 
-3.91(-3) 



-1.61(-3) 
-3.45(-3) 
-6.60(-3) 
-1.02(-2) 
-9.66(-3) 
-3.80(-3) 



-1.5(-4) 

-5.8(-4) 

-1.73(-3) 

-4.01(-3) 

-5.53(-3) 

-2.97(-3) 



-2.7(-5) 

-2.7(-4) 

-1.12(-3) 

-2.53(-3) 

-1.98(-3) 



6.4. Hanle effect with source function and opacity fluctuations 

In the presence of a magnetic field, the problem loses all its symmetries and all twelve unknown amplitudes Ai^S'q 
must be considered. Yet, there are three important cases in which the field is symmetric enough to preserve some 
symmetries, which lead to simplifications of the description and are interesting for understanding the radiative transfer 
problem. We will study in some detail the following three configurations: vertical magnetic field, and horizontal 
magnetic field aligned with, or transversal to, the thermodynamic fluctuations. The simplifications to which these 
symmetric configurations lead are summarized in Table |U 

It is important to note that Ai^Sq — > (for k ^ 0) when approaching to the stellar surface, regardless of the 
magnetic field configuration. Consequently, Ai^SV, Ai^Sq, Ai^Sj/ — > 0, and there are no observable multidimensional 
effects in the limit of tangential observation. In the following, we shall consider the case of finite fj, > observations. 

As is well known, a vertical field produces no Hanle effect in a plane-parallel medium (just consider sm0B = 
and J\ = Jf = in Eqs. ©)• This is no longer true in the presence of horizontal inhomogeneities, since the vertical 
direction is no longer a symmetry axis for the radiation field. Symmetry considerations and Eqs. ([9J for the amplitudes 
Ai,2 Sq show that the only non vanishing amplitudes in this case are the ones listed in Table (j4|). Moreover, when the 
magnetic field is strong enough to reach saturation (r — > oo), the only non vanishing amplitude is Ai5q. Figure li"2l 
shows AQ/I and AU /I emergent at x — (where the fluctuation of the Planck function is maximum), for a ray along 
the y-axis and /i = 0.1 (close-to-the-limb observation), or with fi = 1 (disk center observation or forward scattering 
geometry) . Upper panels show that the effect of a vertical magnetic field is a depolarization of the signal produced 
by the horizontal fluctuations, and a rotation of the polarization plane. For strong fields (r 3> 1), depolarization is 
partial in the fi — 0.1 geometry, and complete at fi = 1. 

When the magnetic field is horizontal and aligned with the LOS (central panels), the effect when observing at 
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Fig. 9. — Center-to-limb variations of amplitudes Ail/I, AiQ/I, and A2U/I corresponding to the amplitudes in Fig. (5), for observation 
centered and along the slabs, i.e., at x = 0, and along the y axis (LOS with \ = 0). Labels indicate wavenumber k. 




Fig. 10. — Emergent fractional polarization Q/I close to the limb at fi = 0.1 (left), and at disk center (right) from an atmosphere with a 
fluctuating Planck function B = B[l + 0.2 cos(fea;)] ; labels indicate wavenumber k. The LOS in the left panel is taken at x = and along 
the slabs (i.e., along the y-axis, see Fig.[TJ, and the positive direction of Q is perpendicular to the limb. The line of sight in the right panel 
is taken at x = and the positive-Q direction is perpendicular to the slab. 

H = 0.1 can be interpreted as in the classical picture in a planeparallel atmosphere (a depolarization, which may 
be total, and rotation of the polarization plane). In forward scattering geometry, the field destroys the orthogonal 
polarization generated by the horizontal inhomogeneities, while no rotation of the polarization plane is possible for 
symmetry reasons. 

When the magnetic field is perpendicular to the LOS, no rotation of the polarization plane is possible at any fi, due 
to symmetry. The magnetic field lies along the polarization direction which leads to partial depolarization at fi = 0.1 
due to the precession of the LOS component of the dipole. In forward scattering (fi = 1), however, overpolarization is 
possible because the anisotropy is negative at t « 1. 

All the previous discusions apply to spatially resolved observations. It is however interesting to consider the effect of 
finite spatial resolution by averaging the fluctuations over intervals A along the x axis. Thus, the average amplitudes 

at x — are (61) — J A ^ 2 81 Ax, with analogous expressions for (SQ) and (SU). The emergent fractional polarization 
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x=P/4 



-2 -n.fi I I -0.1 1 -0.1 
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Av/Av D Av/Av D Av/Av D Av/Av D 

Fig. 11. — Emergent fractional polarization from an atmosphere with horizontal fluctuations of the Planck function and opacity. The 
wavenumber of the fluctuations is k = 0.3, the Planck function amplitude is AB/B = 0.2, and the opacity fluctuations are a = +0.1 (solid 
lines labeled with a + sign), —0.1 (solid lines labeled with a — sign), or (dotted lines). Leftmost panels correspond to a close-to-the-limb 
(fi = 0.1) LOS; rightmost panels to a disc center (/i = 1) geometry. Q/I is shown at x = with a LOS along the ?/-axis; U/I is shown at 
x = P/4 (P = 2n/k is the period of the horizontal perturbation). 



TABLE 4 

Independent unknowns according to the magnetic field geometry 
""Geometry 1 non-vanishing Ai^Sq amplitudes^ 

No magnetic field (B x = B y = B z = 0) A i s ?v A 2 S i> AiSj 

Vertical field (B x = B y =0) AiSg, [A 2 Sj], [A 2 S|], [AiSf], [AiSf] 

Horizontal field across slabs (B y = B z = 0) AiS 1 ^, [A 2 5f], [AiSf], AiSj, [A 2 Sf] 

Horizontal field along the slabs (B x = B z = 0) AiSff, [AgSg], [Ai^Sf], Aigf, [A 2 S|] 

"Horizontal fluctuations ~ cos(kx) (invariant direction along y-axis). 
'Elements between brackets vanish when the field saturates (r — > oo). 

with low spatial resolution at x = is then 

(Q) = Q(*m) + AiQ(^m) sinc(^) <£/) = ^(z M ) + A 1 L/(z M )sinc(^) 

(/> J(2 M ) + Ai/(« M ) sinc(^) ' (/) I{z M ) + *iI(z M ) sinc(^) ' 

where zm is the outer free boundary of the atmosphere and sinc(x) = sin(x)/x. At high spatial resolution, A — > and 
we recover previous results; at low spatial resolution, A — > oo and we recover the fractional polarization Q/I and U/I 
of the average model. Actually, for fcA/2 > n, the amplitudes fall below 20% the resolved value, and for fcA/2 > 47r the 
amplitudes are always below 10% the resolved value. Therefore, for very low spatial resolutions the effects of horizontal 
fluctuations disappear and we recover the plane-parallel limit. In particular, if we have a large-scale magnetic field in 
an atmosphere with horizontal fluctuations of the thermodynamic parameters, in the limit of low resolution only the 
Hanle effect of the large scale magnetic field remains. 

7. CONCLUSIONS 

In this paper we have formulated and solved the problem of scattering polarization in horizontally inhomogeneous 
stellar atmospheres for continuum and resonance line radiation including the Hanle effect from random and determin- 
istic magnetic fields. If the amplitudes of the horizontal fluctuations of the Planck function and of the opacity are 
small enough (compared with the spatially-averaged values), the ensuing radiative transfer problem can be linearized; 
furthermore, if the opacity does not fluctuate, the response of the radiation field to Planck function fluctuations is 
already linear. In both cases, two- and three-dimensional scattering polarization problems can be solved at the compu- 
tational cost of one-dimensional problems, because no spatial discretization is needed along the horizontal directions, 
thus avoiding computationally costly interpolations between grid points. We have developed and implemented iterative 
numerical methods and formal solvers of the transfer equations to solve this type of multidimensional radiative transfer 
problems. 

The numerical codes developed implementing these methods allow carrying out several interesting investigations. 
For instance, it is important to study the impact of the symmetry breaking effects caused by the presence of horizontal 
atmospheric inhomogeneities on the emergent scattering polarization and its importance relative to the modification of 
the linear polarization produced by the Hanle effect. Here we have given a first step towards such a goal by considering 
sinusoidal horizontal inhomogeneities in the Planck function and in the opacity, in the absence and in the presence 
of deterministic magnetic fields (i.e., with a fixed orientation). The most important results of our study, to be taken 
into account when interpreting high-spatial resolution observations of the linear polarization produced by scattering 
processes in the solar atmosphere, are the following: 

(1) When considering increasingly smaller horizontal atmospheric inhomogeneities (i.e., increasingly larger wavenum- 
bers) we find that the amplitudes of the fractional scattering polarization signals first increase and then, beyond a given 
wavenumber, the inhomogeneities become optically thin and the planeparallel limit is recovered again (see Fig. 6). The 
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Fig. 12. — Amplitudes of the emergent Stokes parameters AQ/I and AU/I at x = 0, observed along the x-axis at fi = 0.1 (left panels), 
and at fi = 1 (right panels). A periodic horizontal fluctuation (~ cos(0.3a;)) is assumed for the Planck function. Top panels correspond 
to a vertical magnetic field; middle panels to a horizontal field along the y-axis; bottom panels to a horizontal field along the rr-axis. For 
a given magnetic field configuration the magnetic field strength is T = (solid line), T = 0.3 (long-dashed), T = 1 (dashed), and T = 10 
(dotted). The Q > direction is parallel to the horizon on the left panels and transversal to the "slab" (x-axis) on the right ones. 

maximum impact occurs for /c«0.2 or fc«0.3, which corresponds to P-values of about 20 opacity scale heights. Note that 
the symmetry breaking effects caused by the horizontal atmospheric inhomogeneities may produce forward-scattering 
Q/I and U/I signals, without the need of an inclined magnetic field. 

(2) In the presence of atmospheric inhomogeneities with a dominant horizontal scale of variation (e.g., such as those 
suggested by the fibrils seen in Ha high-resolution images) even a vertical magnetic held may modify significantly the 
linear polarization corresponding to the zero-held reference case. We find that the presence of a vertical field reduces 
the amplitude of the horizontal variation of the emergent Q/I profile and produces U /I signals. This occurs in both, 
the scattering geometry of a close to the limb observation and in the forward scattering geometry of a disk center 
observation. However, in this latter case a vertical field with strength in the Hanle saturation regime kills completely 
the horizontal fluctuation of the Q/I signal, in addition to that of U/I. 

(3) If the magnetic field in the above-mentioned scenario is significantly inclined (e.g., horizontal) and perpendicular 
to the orientation of the "fibrils" we may have an increase in the amplitude of the forward-scattering signals produced 
by the symmetry-breaking effects. 

Finally, we point out that the results given in this paper in graphical and tabular form, and the general symmetry 
and dimensional constrains derived in Sect. l2.4l and l4.3l serve as benchmarks for general radiative transfer codes. With 
ours we plan to investigate how the presence of inhomogeneities in the solar chromosphere may affect interpretations 
in terms of the Hanle effect of the linear polarization signals in spectral lines like Ca II K, Mg II k and Lycv. 

Financial support by the Spanish Ministry of Science and Innovation through projects AYA20I0-I8029 (Solar 
Magnetism and Astrophysical Spectropolarimetry) and CONSOLIDER INGENIO CSD2009-00038 (Molecular As- 
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APPENDIX 

A a/3 OPERATORS AND THEIR SYMMETRY PROPERTIES 

This appendix derives general expressions for the A Qj g operators introduced in equations (|47[) - (|49[) that give the 
explicit dependence of the radiation field tensor components Jq on the statistical elements Sq . There is also a 
discussion on the symmetry properties of these operators according to the dimensionality of the problem. 

The spherical components of the radiation field tensor can be written in a condensed manner as 

f°° f rlQ . 3 

Jq=J^ ^duj—Yj^^mi^), (ai) 

where Ii(v, ft) (i = 0, 1, 2, 3) stands for the Stokes parameters I„n, Qvn, U„n and V u n, respectively, while Tq are the 
polarization tensors introduced by Landi Degl'Innocenti (1984). Analogously, the line source functions for each Stokes 
parameter can be expressed as 

KQ 

where wf^j = 1 and w Jy 2 ^ has been introduced after Eqs. @ (atomic orientation is irrelevant in our problem and 

w^pj is not necessary here). The total source function including the effect of an unpolarized background continuum 
is then 

Si = r„Sf nc + (1 - r v )S cmt S ifi , (A3) 

where 5i, is the Kronecker's delta and r v has been defined after Eqs. ((T|). 
From Eqs. (flj, the Stokes parameters at a point x in the atmosphere can be expressed as follows (e.g., Mihalas 1978) 

[/<(!/, n)].= / S i (x')e^ m >^(x')ds + lP(x t n)e^ m > Wa ) =A va [Si] +li(y,tt). (A4) 
Ji{x,ci) 

where the integral must be understood along the ray path x' = x'(s) with direction fi, from the boundary point Xq 

to x. D ' i s the incident radiation field at Xq and t{x,x') the optical depth between x and x' . The second identity 
in Eq. ()A4[) introduces the integral operator A^o and the attenuated Stokes parameters 2,-, 
From Eqs. (|AT|) - (|A3|) we get 

3 

j« = f Av <t> v I ^ Y, w j K j } E t q (*> ^ r Q ; & °) A ^ ^ 

J 11 K'Q' i=0 

3 

where linearity of the A„f2 operator has been applied. 
Introducing the following A-like operators and attenuated radiation field tensor: 

3 

**q,k>q>[S$] = [du^ <f ^ W ^J2 T Q^ n)TQ>\h n)A„ n [r„ S^}, (A6) 
A KQ AS cont } = Jdv^ j>^-Tg(Q,n)Kn[{l -r„)S cont ], (A7) 

3 

J$= [ cf> u du I^YTQ^ 1 ^^ ( A8 ) 
the radiation field tensor components can be expressed as 

2 K 

J Q = E E ^K Q ,K>Q>lS«;}+A KQ AS cont }+Jg- (A9) 

K'=-2 Q'=-K 

Expressions such as these have been derived by Landi Degl'Innocenti, Bommier & Sahal-Brechot (1990). However, it 
is more convenient from a numerical point of view to obtain similar expressions for the real (and linearly independent) 
components introduced through Eqs. ([2])-([3]). The expressions thus obtained (Manso Sainz 2002), can also be derived 
within the equivalent phase matrix formalism (e.g., Anusha & Nagendra 2011). 
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TABLE 5 

Analytical Expressions for vj a p(Cl) 

rooo = 1 

ro i = roio = ^= (3/i 2 - 1) 
roo2 = 2OT20 = — v3jn ^/ 1 — /j 2 cos x 
TO03 = 2ro30 = v^mV 1 ^ sin X 
ro 04 = 2ro 40 = -^(1 - M 2 ) cos 2x 
ro 05 = 2ro 50 = -^(1 - M 2 )sin2x 



roil 


= \[C^ 2 - 


1) 2 +9(1- M 2 ) 2 ] 




roi2 


= — 2ro2i = 


yf/iV 1 - M 2 cos X[l - 3/i 2 + 


3(1 -M 2 )] 


roi3 


= — 2ro3i = 


yiWl-V sinx[(3 M 2 -l) 


+ 3(^-1)] 


roi4 


= 2ro4i = a 


/|(1 + 3m 2 )(1-M 2 )cos2 X 




tX>15 


= 2ro5i = - 


^y§(l + 3M 2 )(l-M 2 )sin2x 





ro 22 = -§(!- M 2 )(2M 2 cos 2 X + sin 2 x) 

^23 = ^32 = |(2/i 2 - 1)(1 - M 2 )cosxsinx 

ro24 = — ro42 = — f MVl — /i 2 (cosxcos 2x/i 2 + sinxsin2x) 

ro25 = — ro 5 2 = §Mvl — M 2 (cos x sin 2x/i 2 -sinxcos2x) 



W33 = - M 2 )(c os 2 X + 2/i 2 sin 2 X ) 

ro 34 = -ro 43 = f mV 1 ~ M 2 ^ 2 cos2xsinx - cosxsin2x) 

ro35 = — ro 5 3 = — §/4\/i - ^ 2 (^ 2 sin2xsinx + cosxcos2x) 

ro 44 = | [cos 2 2 X (1 + fi 4 ) + 2m 2 sin 2 2 X ] 
ro 4 5 = ro 54 = -|cos2xsin2x(/i 2 - l) 2 

ro 55 = | [2fi 2 cos 2 2 X + (1 + A* 4 ) sin 2 2 X ] 



To this end, we write the first term on the right hand side of Eq. (|A5|) 

3 

j Av(f>v /^X;7^(i,n)A^[r„ 'EwffiTffrtyS?' +rv ]T E 4^ fi)^' +T5 Q ,(i, 

^ i=0 AT' if' g'>o 

(AfO) 

Since the Tq(i,£l) tensor satisfies the conjugation property T^q = {~~^)®\Tq\*i then T^q/ (i, J1)55q/ = 
[7qT («, J^S'q/ ]*■ Let 7q and 7q be the real and imaginary parts of Tq , respectively. The last bracket in equation 
(|A10|) then reads 

T$'(i,tt)s$ +T K Q ,(i,n)s K Q , = 2Re{T£\i,n)s$} = 2 [f«;'(i,n)s$ - T$\i,n)s$]. (Aii) 

Finally, using this expression in Eq. (|A10|) . and introducing the symbol cq> = 2—5q'q, we obtain the explicit dependence 
on Sq and Sq of the real and imaginary components Jq and Jq of the radiation field tensor: 

3 3 
J J K'Q'>Q i=0 i=0 (A12) 

+ J to fa j ^fef (0)A^[(1 - r„) S cont ] + J«, 

3 3 

K'Q'>0 t=0 l= o (A13) 

J to fa j ^7ef (0)A„ n [(l - r„) S cont ] + j§, 
which, enlightening the notation, can be written as 

Jq = E k kq,vq^ + E A ^,^[^ ; ] + Vg.r'i + ^o. ( A14 ) 

K'Q' K'Q' 
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TABLE 6 



Analytical Expressions 
for ro ac (fl) 



roOc = 




roi c = 




OT2 C = 


— sin 9fi cos x 


OT3c = 


— sinS^tsinx 


TO4 C = 


^(l- A1 2 )cos2 X 


rose = 


^(l- M 2 )sin2 X 



Aoo 


Ao, 


A 02 


A 03 


A 4 


A 05 


A, 


A„ 


A 12 


A, 3 


A 14 


A 15 


A 20 


A 2 , 


A 22 


A 23 


A 24 


A 25 


A 30 


A 31 


A 32 


A 33 


A34 


A35 


A 40 


A 4 , 


A 42 


A43 


A 44 


A 45 


A 50 


As, 


A 52 


A 53 


A 54 


A55 



Fig. 13. — Non-vanishing operators according to the medium's geometry. Dark blocks: plane-parallel case. Grey blocks: Cartesian 
two-dimensional medium. In a general three-dimensional atmosphere all of them are non-zero in general. 



E A KQ,KV'^Q'^ + E K KQ.K'Q'y u Q 
K'Q' K'Q' 



8$] 



where 



A 



KQ,K> 
A 



KQ,c 



dv(j)u <p— CQI w~ ~, Knr v , 
<l,/ "<- $ ~j~T w kq,c A*/n(l - r v ) , 



(A15) 

(A16) 
(A17) 



and analogously for the remaining operators. The angular weights are 



7q (i)T£'(i), 



KQ.K'Q' 



KQ.K'Q 1 



3 

E 

i=0 
i=0 



i=0 
3 



^,0 = ^(0). 

^ = 7^(0). 



Their explicit expressions are given in Tables [5] and [5] after a convenient index renaming (0, 1, 2, 3, 4 and 5 for 00, 20, 



and 22, respectively). Thus, we get (cf. Eqs. (|47| 


)■ 












= A 00 SS° - 


hA 01 5S^- 




h A 03 <5^ + 


A 04 <^ - 


h A 05 8S% 


+ A ac 8S cont 


+ SJl 


(A18a) 


6J* 


= Aio^^ " 


h AuSSq - 


hA 12 W 2 - 


h A 13 5S? + 


A 14 8S% - 


h A 15 5S* 


+ A lc 8S cont 


+ SJl 


(A18b) 


8 J* 


= A 2 o8S° - 


hA 21 5S£- 


hA 22 W 2 - 


h A 23 5^ + 


A 2A 8Sl - 


h A 25 5S* 


+ A 2c 8S cont 


+ SJ^ 


(A18c) 


8Jl 


= A 30 8S° - 


bA 31 8S%- 




h A 33 5S? + 


A 34 8S* - 


h A 35 5S| 


+ A 3c SS cont 




(A18d) 


SJ% 


= A^qSSq - 


hA 41( 55 2 - 


hA 42 W 2 - 


h A 43 <5S? + 


A448S2 - 


- A^SS^ 


+ A Ac 8S cont 


+ SJl 


(A18e) 


8Jl 


= A 50 5S% - 




^A 52 8Sl- 


h A 63 5^ + A 54 5^ - 


h A 55 <5^ 


+ A^SS^ 




(A18f) 



where the matrices A a p are obtained from A n (see Eqs. (|48|) -(l49 
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In a plane-parallel atmosphere, the A„n(i, j) elements are independent of x an d the azimuthal integrals in Eqs. (| A16[l - 
(IA17J) can be performed analytically. It is then found that all the non-diagonal operators but Aoi and Am, vanish, as 
well as A QC (a > 2). Furthermore, A22 = A33 and A44 = A55. 

In a two-dimensional atmosphere, taking the y-axis as the direction invariant under translations, the A„n operator 
satisfies the symmetry relation 

A w ,( M , x )(i,j) = A Ut(tfl - x )(i,j). (A19) 

Taking into account relation (|A19|) in Eqs. (IA16|) - (|A17|) . we find that A 03 , A 5, A 13 , A 15 , A23, A25, A 34 , A 45 and 
their symmetric ones vanish, as well as A 3c and As c . 

These results are summarized in Fig. 1131 Dark blocks show the non-zero operators in a plane-parallel atmosphere; 
in a two-dimensional medium gray blocks are non-zero too. Note that the Sf and Sf are radiatively coupled only 
between them. Therefore, in the absence of magnetic field couplings, they are decoupled from the thermal source and 
vanish everywhere in the atmosphere. This is another proof of the result obtained in §4.3 from general symmetry 
considerations. In a three-dimensional medium, all the operators are, in general, non-zero. 
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